import cmocean as cmo
import numpy as np
import scipy.constants as spc
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation
from IPython.display import display, Math, HTML
from matplotlib.animation import FuncAnimation
plt.rcParams.update(
{
'mathtext.fontset':'cm',
'font.family':'serif',
'font.sans-serif':'Times New Roman',
'figure.figsize':[9,4],
'figure.titlesize':22,
"figure.dpi":90,
'savefig.dpi':90,
'axes.titlesize':20,
'axes.labelsize':18,
'axes.titley': 1.0,
'axes.titlepad': 5.0,
'axes.edgecolor':'black',
'axes.grid': False,
'grid.alpha': .5,
'xtick.labelsize':14,'ytick.labelsize':14,
'xtick.major.size':6,'ytick.major.size':6,
'xtick.major.width':1.25, 'ytick.major.width':1.25,
'xtick.direction':'inout','ytick.direction':'inout',
'xtick.top':False, 'ytick.right':False,
'legend.title_fontsize':14, 'legend.fontsize':14,
'legend.borderaxespad': 1, 'legend.borderpad': 0.5,
'legend.framealpha': 1,
'legend.handleheight': 0.5, 'legend.handlelength': 2.0, 'legend.handletextpad':0.5,
'legend.labelspacing': 0.25,
'legend.fancybox':False,
'legend.edgecolor': '0',
'legend.frameon': True,
'legend.markerscale': 1.25,
'animation.embed_limit':2**128,
'animation.html': 'jshtml'
}
)
pi = spc.pi
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
default_plots()
def periodic(f):
f[0,:] = f[-2,:]
f[-1,:] = f[2,:]
f[:,0] = f[:,-2]
f[:,-1] = f[:,2]
return
def fixed(f):
f[0,:] = 0
f[-1,:] = 0
f[:,0] = 0
f[:,-1] = 0
return
c = 1
x = np.linspace(0,1,128)
y = np.linspace(0,1,128)
xx, yy = np.meshgrid(x,y)
u0 = np.zeros(y.shape+x.shape)
steps = 500
dx = x[1]-x[0]
dy = y[1]-y[0]
dt = (dx/c)/2
t = np.linspace(0,steps*dt,steps+1)
u = np.zeros(t.shape + u0.shape)
u_t = np.zeros(t.shape + u0.shape)
u[0] = u0/10
u_t[0] = np.zeros_like(u0)
rx = (c/dx)**2
ry = (c/dy)**2
# def source(f, n_t):
# return np.cos(f*pi*n_t/steps)**2
xL = .5-dx
xR = .5+dx
yB = .5-3*dy
yT = .5+3*dy
mask = ~(((xx < xR)*(xx > xL)*(yy>yT)) + ((xx < xR)*(xx > xL)*(yy<yB)))
source = (xx-.4)**2 + (yy-.5)**2 < .002
boundary = 'absorbing'
src_int = 40
stp_int = 2
if boundary == 'periodic':
for n in range(steps):
if (n%src_int==0)*(n < steps//stp_int):
u[n] += source
u = u*mask
periodic(u[n])
u_t[n+1,1:-1,1:-1] = u_t[n,1:-1,1:-1] \
+ dt*rx * (u[n,2:,1:-1] + u[n,:-2,1:-1] - 2*u[n,1:-1,1:-1]) \
+ dt*ry * (u[n,1:-1,2:] + u[n,1:-1,:-2] - 2*u[n,1:-1,1:-1])
u[n+1,1:-1,1:-1] = u[n,1:-1,1:-1] + dt*u_t[n+1,1:-1,1:-1]
if boundary == 'fixed':
for n in range(steps):
if (n%src_int==0)*(n < steps//stp_int):
u[n] += source
u = u*mask
fixed(u[n])
u_t[n+1,1:-1,1:-1] = u_t[n,1:-1,1:-1] \
+ dt*rx * (u[n,2:,1:-1] + u[n,:-2,1:-1] - 2*u[n,1:-1,1:-1]) \
+ dt*ry * (u[n,1:-1,2:] + u[n,1:-1,:-2] - 2*u[n,1:-1,1:-1])
u[n+1,1:-1,1:-1] = u[n,1:-1,1:-1] + dt*u_t[n+1,1:-1,1:-1]
if boundary == 'absorbing':
for n in range(steps):
if (n%src_int==0)*(n < steps//stp_int):
u[n] += source
u = u*mask
u[n+1,1:-1,0] = (1-c*dt/dx)*u[n,1:-1,0] + (c*dt/dx)*u[n,1:-1,1]
u[n+1,1:-1,-1] = (1-c*dt/dx)*u[n,1:-1,-1] + (c*dt/dx)*u[n,1:-1,-2]
u[n+1,0,1:-1] = (1-c*dt/dy)*u[n,0,1:-1] + (c*dt/dy)*u[n,1,1:-1]
u[n+1,-1,1:-1] = (1-c*dt/dy)*u[n,-1,1:-1] + (c*dt/dy)*u[n,-2,1:-1]
u_t[n+1,1:-1,1:-1] = u_t[n,1:-1,1:-1] \
+ dt*rx * (u[n,2:,1:-1] + u[n,:-2,1:-1] - 2*u[n,1:-1,1:-1]) \
+ dt*ry * (u[n,1:-1,2:] + u[n,1:-1,:-2] - 2*u[n,1:-1,1:-1])
u[n+1,1:-1,1:-1] = u[n,1:-1,1:-1] + dt*u_t[n+1,1:-1,1:-1]
def animate2D_mesh(u, t, fps = 24, frames = 96, cmap = cmo.cm.deep):
t_steps = int(len(t)-1)
i_frames = np.arange(0, t_steps, int(t_steps/frames))
fig, ax = plt.subplots(figsize = (7,7), constrained_layout = True)
image = ax.imshow(u[0], cmap = cmap, vmin = np.min(u[0])/2, vmax = np.max(u[0])/2)
ax.axis(False)
def animate(i):
image.set_array(u[i])
plt.close()
anim = FuncAnimation(fig, animate, interval = int(1e3/fps), frames = i_frames)
return anim